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We consider particles suspended in a randomly stirred or turbulent fluid. When effects 
of the inertia of the particles are significant, an initially uniform scatter of particles can 
cluster together. We analyse this 'unmixing' effect by calculating the Lyapunov exponents 
for dense particles suspended in such a random three-dimensional flow, concentrating on 
the limit where the viscous damping rate is small compared to the inverse correlation time 
of the random flow (that is, the regime of large Stokes number). In this limit Lyapunov 
exponents are obtained as a power series in a parameter which is a dimensionless measure of 
the inertia. We report results for the first seven orders. The perturbation series is divergent, 
but we obtain accurate results from a Pade-Borel summation. We deduce that particles 
can cluster onto a fractal set and show that its dimension is in satisfactory agreement with 
previously reported in simulations of turbulent Navier-Stokes flows. We also investigate the 
rate of formation of caustics in the particle flow. 
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1 Introduction 

1.1 Overview 

This paper discusses small particles suspended in a randomly moving fluid. We assume that 
the fluid flow is mixing, and concentrate on cases where the fluid flow is incompressible (al- 
though compressibility is considered too, for completeness). At first sight, it seems as if the 
particles suspended in an incompressible mixing flow should become evenly distributed. If the 
particles were simply advected by the fluid, this is indeed what would happen. However, it has 
been noted that when the effects of finite inertia of the suspended particles are significant, the 
particles can show a tendency to cluster. This remarkable 'unmixing' effect was discussed in a 
theoretical paper by Maxey [T], who proposed that suspended particles (assumed denser than 
the fluid) cluster because they are centrifuged away from vortices. There have been many other 
theoretical papers on this phenomenon, and an experimental demonstration was reported by 
Eaton and Fessler [2]; the literature is reviewed in section [L~2l below. The suspended particles 
are characterised by the rate 7 at which their velocity relative to the fluid is damped due to 
viscous drag, and the random velocity field is characterised by a correlation time r. The product 
= 7T is a dimensionless parameter: in much of the literature St = 1/S1 is termed the 'Stokes 
number'. There is a consensus that the clustering effect is observed when the Stokes number is 
of order unity. 

We investigate a random-flow model with short correlation time, which is susceptible to 
mathematical analysis. We show that in general this model can exhibit pronounced clustering 
in circumstances where St 3> 1, when the centrifuge mechanism is not effective. When the 
model is applied to fully-developed turbulence it predicts that the clustering is strongest when 
St = 0(1) (in agreement with most numerical investigations), but it indicates that the centrifuge 
effect is not essential to understanding the phenomenon. 

Our principal results are series expansions for the Lyapunov exponents of the trajectories of 
the suspended particles in terms of e = K,VSt, where the dimensionless parameter k (defined in 
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section [L3l below) is 0(1) for fully developed turbulent flow, but may be small for other types 
of random flow. Our analysis of the random-flow model is exact in the limit as k — > 0: note that 
in this limit St 3> 1 when e = O(l). We use the Lyapunov exponents to estimate the 'dimension 
deficit' A, which is the difference between the spatial dimension and the (Lyapunov) fractal 
dimension di, of the set onto which the particles cluster: di, = 3 — A (this will be explained in 
section [L~2jh Figures [T]a and[]J> illustrate our results. Figure [l]a compares the value of A as a 
function of e obtained by simulation of our model with results from a Borel summation of our 
perturbation series. The results show that there is fractal clustering when e = 0(1) in the limit 
as k, — ► 0, indicating that fractal clustering can occur for large Stokes numbers. Despite the fact 
that our perturbation series is divergent, we obtain satisfactory results from Borel summation. 
Figure [T)b compares the results from a direct numerical simulations of particles in a turbulent 
Navier-Stokes flow (these data (□) are taken from [3]) with results from our random-flow model: 
the data used for figure [1^ are re-plotted with the value of n chosen to give the best fit between 
the two curves, as judged by eye. (We cannot determine the scaling theoretically because the 
value of k for fully developed turbulence is not known) . The theoretical curve in figure [T]b is 
obtained from our random-flow model, which has vanishingly small correlation time (because 
for given e, St —* oo as k —* 0). The centrifuge effect cannot therefore cause clustering in 
this random-flow model. Our model has a maximum dimension deficit of approximately 0.35, 
as opposed to 0.40 for particles in a Navier-Stokes flow, and the form of the curves is similar. 
We conclude that the random-flow model provides a satisfying degree of agreement with full 
Navier-Stokes dynamics, despite the fact that the centrifuge mechanism plays no role. 

Another mechanism for clustering is the formation of fold caustics in the flow of the particles. 
We show that caustics are prevalent when e > 1. We also consider an explicit expression for 
the sum of the Lyapunov exponents when the Stokes number is small. Taken together with 
earlier results on the overdamped limit (referred to in sections 11.21 El and [9]) , our results give a 
satisfyingly complete understanding of the local dynamics of suspended particles. 

Some of the results were summarised in an earlier letter [4j. 

1.2 Discussion of earlier work 

We first briefly review earlier work on clustering in random flows, before giving a fuller descrip- 
tion of our results. 

Maxey [1] expressed the trajectory of an inertial particle in terms of a 'synthetic' velocity 
field, which is obtained as a perturbation of the velocity field of the fluid. He showed that this 
synthetic velocity field has negative divergence when the vorticity is high or the strain-rate low, 
and predicted that particles would have low concentrations in regions of high vorticity due to 
this 'centrifuge effect'. A correlation between particle density and vorticity has been seen in 
direct numerical simulation of particles suspended in a fully-developed turbulent flow O [6] . 

The mechanism proposed by Maxey led to a prediction concerning the observability of this 
'preferential concentration' effect. It is plausible that the 'centrifuge mechanism' will be less 
effective when the particles are overdamped (fl 3> 1) or when the velocity fluctuates too rapidly 
to allow the density of suspended particles to respond (when <C 1). This leads to the hypothesis 
that the preferential concentration effect should only be observed when f2 is close to unity. 
Experimental work on particles suspended in fully developed turbulent air flows [2] appears to 
support this hypothesis, as do computer simulations [6], [3] . 

An alternative approach arises from work of Sommerer and Ott [8], who discuss patterns 
formed by particles floating on a randomly moving fluid. They characterise these patterns in 
terms of their fractal dimension, and suggest that the fractal dimension can be obtained from 
ratios of Lyapunov exponents Aj of the particle trajectories, using a formula proposed by Kaplan 
and Yorke [9]. The argument of Sommerer and Ott extends directly to particles suspended in 
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turbulent three-dimensional flows. The spatial particle flow is characterised by three Lyapunov 
exponents, Ai > A2 > A3, and provided Ai + A2 > (which is always true for the particles 
moving in an incompressible fluid flow), the Kaplan- Yorke estimate for the fractal dimension is 
determined by the dimensionless quantity 

A = -^(A! + A 2 + A 3 ) (1) 

l A 3| 

which we term the 'dimension deficit'. When A > 0, the Kaplan- Yorke estimate of the dimension, 
termed the Lyapunov dimension, is 

d L = 3 - A (2) 

and c^l = 3 if A < 0. Clustering effects are significant if the fractal dimension is significantly 
lower than the dimension of space. The relations ([I]), ([2]) give a strong motivation to investigate 
the Lyapunov exponents of suspended particles. Bee [TOl [11] has performed detailed numerical 
investigations for a specific ensemble of random flow fields, showing how the Kaplan- Yorke fractal 
dimension varies as a function of the Stokes number for his model flow, reaching a minimum 
for a value of f2 which is of order unity. These numerical calculations of the dimension deficit 
are complementary to the ideas proposed by Maxey in p], in that they quantify the clustering 
effect without explaining its origin. 

The motivation for investigating Lyapunov exponents can also be explained without referring 
to the fractal dimension. For three-dimensional flows, the sum of the three largest Lyapunov 
exponents of the particle trajectories may be negative, implying that volume elements almost 
always contract. However, for incompressible flows the first Lyapunov exponent is always posi- 
tive, implying that nearby particles almost surely separate. If (Ai + A2 + A3)/Ai is negative, there 
is a tendency for particles to cluster, but if the magnitude of this quantity is small compared 
to unity the clustering effect will be negligible, because in that case clusters are stretched and 
folded more rapidly than their density accumulates. 

Most of the theoretical work on clustering in turbulent flows has emphasised instantaneous 
correlations between vortices and particle-density fluctuations. There is an alternative viewpoint 
on the origin of density fluctuations which is also theoretically tenable, and which is the basis for 
our own analysis. It can be argued that the density fluctuations are generated by a multiplicative 
random process: volume elements in the particle flow are randomly compressed or expanded, 
and the ratio of the final density to the initial density after many multiples of the correlation 
time r can be modeled as a product of a large number of random factors. According to his 
picture, the density fluctuations will be a record of the history of the flow, and may bear no 
relation to the instantaneous disposition of vortices when the particle density is measured. The 
particle density is expected to have a log-normal probability distribution, and the mean value 
of the logarithm is related to the Lyapunov exponents of the flow. This is another motivation 
for calculating Lyapunov exponents. 

The idea of considering clustering due to random flows has been used in earlier work which 
adapted results relating to purely advective flows. Most of this literature treats the limiting 
case where the correlation time of the velocity field is very short, a case which is known as 
the Kraichnan model |12j . The Lyapunov exponents of such advective flows have been calcu- 
lated in different ways by several authors: these calculations include results for compressible 
and solenoidal (incompressible) flows: the earliest calculation appears to be by LeJan [13], who 
treated a spatially correlated Brownian motion. Later work [HI [15] subsequently showed that 
the same results apply to a flow with a smooth time dependence, but a very short correlation 
time. These calculations on purely advective flows cannot explain clustering of particles sus- 
pended in incompressible flows, because the density of advected particles remains constant for 
incompressible flow. Elperin [16J proposed to analyse the motion of inertial particles by ap- 
plying results derived for advective flows to the synthetic velocity field derived by Maxey pQ, 
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which has a compressible component. The same approach was subsequently adopted in \17\ . 
The approximations employed in these papers make use of results from the Kraichnan model of 
passive scalars [12], valid in the limit where the correlation time r approaches zero (St — > oo). 
Yet Maxey's synthetic velocity field is obtained in the overdamped limit where, by contrast, the 
Stokes number is small, St — > 0. 

In summary, the clustering effect can be characterised by calculating the Lyapunov exponents 
of the particle trajectories, but there is a limited theoretical understanding of these, based on 
Maxey's correction to advective flow (which is derived in the limit where < 1). There is a 
consensus that significant clustering only occurs when 0, ~ 1 but there is scope for revising this 
expectation. 

Our own work [18\ \19\ [20], [2~Tl 0] uses a model for a random flow with a very short correlation 
time, but the effects of particle inertia are properly accounted for. We remark that a similar 
approach was proposed by Piterbarg [22], who studied the largest Lyapunov exponent in a 
two-dimensional flow. 

1.3 Plan of paper and summary of results 

In section [21 we discuss the dimensionless parameters of the model, and in particular we note 
the significance of an additional dimensionless parameter, the 'Kubo number' (this term is used 
in the plasma-physics literature, [23J). In terms of a typical fluid velocity u and correlation 
length £, the Kubo number is k = ut/£. We argue that k cannot be large, that k = O(l) for 
fully-developed turbulent flows, and that k may be small in other systems, such as randomly 
stirred fluids. 

Most of our paper is concerned with the case where the Stokes number St = 1/Q is large. This 
is the underdamped limit where inertial effects are most likely to be significant and where very 
little work has been done previously. In section [3] we show how the three Lyapunov exponents for 
inertial particles can be obtained from expectation values of a system of nine coupled stochastic 
differential equations. 

In section U] we discuss how these general equations can be represented by a system of 
Langevin equations in the limit where the Stokes number is large and the Kubo number is 
small. A dimensionless parameter e oc kVL~ 1 / 2 plays a natural role in these Langevin equations: 
inertial effects are significant when e is large. 

Section [5] shows how the corresponding Fokker-Planck equation can be mapped to a per- 
turbation of a nine-dimensional isotropic quantum harmonic oscillator. Using the algebra of 
harmonic-oscillator raising and lowering operators, we develop perturbation expansions for the 
Lyapunov exponents to large orders in e, with exact expressions for the coefficients. These 
perturbation expansions are presented in section [6] and are the principal results in this paper. 

In the case where inertia is important it is possible for the density to diverge due to the 
formation of caustics, which are surfaces on which the Jacobian determinant of the particle flow 
field vanishes. In section [6] we also quantify the rate of formation of these caustic surfaces. 
Caustics can have a very significant effect on the aggregation of suspended particles, because 
they can produce a divergent density of particles in a finite time [20 j and because they greatly 
increase the relative velocity of suspended particles [24J . 

The perturbation series are divergent, and in section [6] we also discuss their summation by 
Borel-Pade methods. We find excellent agreement with numerical simulations in some cases, but 
there are revealing discrepancies in other cases: some of our Borel-Pade summations differ from 
Monte-Carlo evaluations by quantities which have a non-analytic dependence on the perturba- 
tion parameter e, of the form exp(— <I>*/e 2 ) (for some constant <$*). In section [7] we consider 
a WKB approach to solving the Fokker-Planck equation, and indicate how such non-analytic 
contributions arise. 
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Section [8] makes connections between our results and earlier work on Lyapunov exponents of 
advected particles: we show that the leading-order term in our perturbation series agrees with 
the Lyapunov exponent for advective flow, and we show that in the limit as k — > the same 
expression for the Lyapunov exponent holds, regardless of whether Q is large or small. In the 
case where the flow is incompressible the sum of the Lyapunov exponents for advected particles 
is zero. In section [9] we calculate the leading contribution due to the sum of the Lyapunov 
exponents for particles in an incompressible flow when f2 3> 1 and k <C 1: the results confirm 
that c^l is very close to 3 when f2 3> 1. 

Our results have important consequences for the theory of particle clustering. The existing 
consensus favours the view that clustering only occurs when ~ 1, and is the result of the 
'centrifuge effect'. Our results present a different picture. Clustering onto a fractal set occurs 
when A, defined by ([I]), is positive, and becomes significant when this number is of order unity. 
Both simulations and the Borel-Pade summations in section [6] indicate that when J] < 1 and 
k -C 1, A is positive when e is less than a critical value e c , and achieves its maximum value of 
A max ~ 0.35 when e equals e max 0.64e c . Thus we establish that clustering can be significant 
when Q <C 1, provided that k<1. Our results on the overdamped limit obtained in section [9] 
indicate that although A > when O, 3> 1, it is very small, implying that clustering effects are 
hard to observe when fi > 1. For fully developed turbulence we have k ~ 1, which is on the 
border of the region of validity of our theory, but a plot of A versus St for our model (figure 
[T]b) shows satisfying agreement with numerical simulations of turbulent flows (after scaling St 
to account for uncertainties in the definitions of correlations times). This indicates that the 
centrifuge mechanism makes a marginal contribution to the clustering process. Systems such 
as randomly stirred fluids, or particles falling under gravity through fully developed turbulence, 
can exhibit flows where fc<l, and may exhibit clustering for large values of St. 



2 Formulation of the problem 

2.1 Equations of motion 

We assume that the particles suspended in the fluid flow satisfy the equations of motion 

r = — p , p = mj[u(r, t) — f] (3) 
m 

where r = (ri,r2,rs) is the position of a particle, p is the particle momentum, m is its mass 
and u(r, t) denotes the velocity field. We neglect effects due to the inertia of the displaced 
fluid: this is justified when p p /pf 3> 1, where pf, p p are the densities of the fluid and particles 
respectively. Equations ([3]) are appropriate for spherical particles when the Reynolds number of 
the flow referred to the particle diameter is small, and when the particle radius a satisfies a <C £. 
Further conditions are required for the validity of this formula: these can always be satisfied if 
the radius of the particle and the molecular mean free path of the fluid are sufficiently small 

Stokes's formula gives the relaxation rate 

Girapfis 9pfU 

7 = = o 2 4 

where v is the kinematic viscosity of the fluid. The neglect of the mass of the displaced fluid 
is an excellent approximation for aerosol systems, and satisfactory for many examples of solid 
particles in water. 

We also assume that the effect of Brownian diffusion of the particles is negligible: the ratio 
of the particle diffusion constant to the molecular diffusion constant of the fluid is proportional 
to the ratio of the molecular mean free path to the particle diameter. 
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2.2 Dimensionless parameters 

We characterise the random velocity field by its statistics, denoting the expectation value of a 
quantity X by {X). We assume that the mean velocity is zero: (u(r, t)) = 0. The velocity 
field can be characterised by its correlation function, which has a correlation length £ and a 
correlation time r. The fluctuations of the fluid velocity have a characteristic scale u: in a 
single-scale flow we would define u 2 = (u 2 ), but in fully developed turbulence it is more natural 
to define u as the velocity scale associated with fluctuations on the dissipative scale, u = (fz^) 1 / 4 , 
where £ is the rate of dissipation per unit mass [26J. 

The equations of motion ([3]) are characterised by four dimensional parameters. The random 
velocity field is described by three scales: u, £ and r. In addition, the interaction of the fluid 
with the particles is described by the damping rate 7. (The mass m can be eliminated from 
the two components of ([3]), but may appear in expressions which contain forces). From the four 
quantities £, r, u, 7 we can form two independent dimensionless groups, the Kubo number, k and 
the Stokes number, I/O. The degree to which the velocity field is compressible is described by 
a further dimensionless variable, T, which will be defined in section HI The average number of 
suspended particles per unit volume, No, is associated with a further independent dimensionless 
parameter, T = Nq£ 3 . The set of dimensionless parameters of the system is therefore 

« = y, ^ = 7T, T, T = iVo£ 3 . (5) 

We consider flow fields ranging from incompressible flow (which corresponds to T = 2, see 
equation in section H]) to pure potential flow (r = g). We concentrate on the underdamped 
limit (fi< 1, large Stokes number), but also give results for the overdamped situation. 

We argue on physical grounds that k cannot be large, and that k ~ 1 for fully developed tur- 
bulence. For real flows satisfying the Navier-Stokes equations, the velocity field is self-convected, 
so that temporal variations of the velocity field at any point are partly due to the 'sweeping' 
action of the flow. If u is the characteristic velocity and £, r are the spatial and temporal corre- 
lation scales, then the transport of the velocity field by its own action will cause it to fluctuate 
on a time scale £/u, which cannot be less than the actual correlation time of the field. Thus 
k = ut/£ cannot be large. Small Kubo numbers are realised in randomly stirred fluids where 
the Reynolds number is small enough that the flow does not spontaneously generate turbulence. 
For fully developed turbulence, the Kolmogorov theory [26j indicates that u, £ and r are all 
functions of £ and u, and dimensional considerations imply that k ~ 1. Our analytical results 
are all derived in the limit where k <C 1, so that fully-developed turbulence is on the borderline 
of applicability of our theory. 

A practically important measure of the degree of clustering is given by the ratio of the 
largest observed particle density -/V max to the mean particle density. The parameter T will have 
a pronounced effect on the distribution pattern of the particles. When T is sufficiently small, 
the particles will appear as a random scatter if the largest Lyapunov exponent (defined below) is 
positive. The set of particle positions is a point set which randomly samples a fractal measure, 
but the fractal is only visible when T is sufficiently large. For large T, the density-enhancement 
factor N max /No may be very large, even when the parameter A (defined by (pQ)) is small. A 
quantitative discussion of these issues would be quite lengthy. Accordingly, in this paper we 
will consider only locally defined properties of the particle trajectories, namely the Lyapunov 
exponents and the rate of caustic formation, both of which are defined below. 

We remark that the dimensionless parameters k and T have apparently not been considered 
in earlier papers on clustering of particles in random flow fields. Uncertainties about the intended 
values of these parameters makes some of the literature quite difficult to understand. 
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2.3 Definitions of rate constants 



The Lyapunov exponents Aj, i = 1, 2, 3 are rate constants which are defined in terms of the time 
dependence of small separations of trajectories from a reference trajectory r(t). We consider 
three trajectories which have infinitesimal displacements from the reference trajectory, <5r,(i), 
i = 1,2, 3. We then consider the length Sr = \5r\\ of a small separation between two trajectories, 
the area 5 A = \r% A r%\ of a parallelogram spanned by two separation vectors and the volume 
5V = Irv^ArsI of a parallelepiped spanned by a triad of separations. The Lyapunov exponents 
are defined by writing 



If Ai < 0, pairs of particles coalesce with probability unity If Ai > and Ai + A2 < 0, particles 
cluster onto randomly moving lines, which stretch and fold. If Ai + A2 > but Ai + A2 + A3 < 0, 
the particles cluster onto randomly stretching and folding surfaces. 

2.4 Caustics 

Another locally defined statistic is the rate of caustic formation. There is no constraint which 
prevents two of the three vectors defining the separation between nearby particles becoming 
collinear, so that the volume element 5V becomes zero for an instant in time. These events 
correspond to 'caustics', where faster moving particles overtake slower ones |27[ l20l [2^] (see 
figure [2] for an illustration in one spatial dimension). Caustics influence the spatial particle 
distribution and the relative velocities of nearby particles. 

There is an increased density of particles on the fold caustics (which are a pair of points 
in the one-dimensional example of figure El but which form surfaces in three dimensions) , and 
the particle density on the caustics diverges in the limit as T — ► 00. This effect is discussed in 
[271 [20] : it is analogous to the divergence of light intensity on optical caustics [28] . 

The other effect of caustics is that the particle velocity field becomes multi-valued in the 
region between the caustics (in figure [2] the velocity field is triple- valued between the caustics). 
Because particles at the same position are moving with differing velocities, the caustics enhance 
the rate of collision of suspended particles |27[ [20] (this has no analogue in optical caustics). 
Caustics therefore facilitate the aggregation of suspended particles. 

We define J, the rate of caustic formation, as the rate at which events where 5V = occur 
for a given triplet of nearby trajectories. We define a dimensionless rate J' by J' = J/7. 

3 Equations determining the Lyapunov exponents 

3.1 Stochastic differential equations for the Lyapunov exponents 

Linearising the equations of motion ([3]) gives 



Ai + A2 + A3 



Ai + A 2 



Ai 



lim - log e ((5r) 

t— >oo t 

lim -log e (SA) 

f^oo t 

lim -log e (JV) . 

t— >oo t 



(6) 



Sp = —^5p + F(t)6r 

5r = —5p . 

m 

where F(t) is a matrix with elements proportional to the rate-of-strain matrix: 



(7) 




(8) 
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(from this point we will use Greek subscripts to label components in three dimensional space, 
reserving Roman indices for components in a nine-dimensional space which appears later). To 
determine the Lyapunov exponents we consider three trajectories displaced relative to a reference 
trajectory by (5r^,5p^), with fj, = 1,2,3. We choose to parametrise the spatial displacements 
as follows 

Sri = -^ini 

5r 2 = X 2 (ni + 59x12) 

5r 3 = X 3 [m + 595<P(Zn 2 + n 3 )] (9) 

where the n^(t) form a triplet of orthogonal unit unit vectors. Nine variables are required to 
parametrise the spatial displacements Sru,- Two parameters specify the direction of ni, and a 
further angular parameter specifies the direction of n 2 relative to ni. There is only a binary 
choice in the direction of ri3, and we resolve this by requiring continuity. This means that a 
further six parameters are required, and ([9]) does indeed contain a further six parameters, namely 
(X 1 ,X 2 ,X 3 ,Z,69,5<j>). 

For the choice of parametrisation above, the small elements required in ([6]) are 5r = X\, 
5A = X\X 2 59 and 5V = X\X 2 X 3 59 2 5(j). We can then extract the Lyapunov exponents from 
expectation values of logarithmic derivatives of the variables used in our parametrisation: 



A 



Xi 
Xi 



X 3 / \5<t>/ \59 

The angles 59 and 5<p both decrease with probability unity, and eventually we retain only leading 
orders in these variables. Initially, however, we retain all terms. 

The equations of motion for the displacements 5r^ also contain corresponding increments of 
momentum, 5p^. These will be parametrised in terms of the nine elements of a 3 x 3 matrix R 
by writing 

S Pfl = KSr^. (11) 
Substituting this relation into ((ZJ) gives the equation of motion for R: 

^ = -7R"-R 2 + F(t). (12) 
at m 

Now we determine the equations of motion for the parameters determining 5r. Differentiating 
each line of © with respect to time, and then using the second equation of ([7]), we obtain 

Xini + Xihi 
1 

= -JiR t ni 

m 

X 2 (ni + 59n 2 ) + X 2 (ni + 59h 2 ) + X 2 59n 2 

= — X 2 R(t)(m +59n 2 ) 
m 

X 3 [m + 595<p{Zn 2 + n 3 )] + X 3 [ht + 695<p(Zn 2 + ri 3 )] 
+X 3 (595(j) + 595<j))(Zn 2 + n 3 ) + X 3 595<pZn 2 

= -X 3 K(t)[ni + 595<P(Zn 2 + n 3 )] . (13) 
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We now take the scalar product of each of these three equations with each of the unit vectors 



n M in turn. It is convenient to use the notations 

R'^(t) = n„(i) • R(t) n„(i) , F'^(t) = n M (t) • F(t) n,(i) (14) 

for the elements of the tensors F and R transformed to the rotated basis. Taking the scalar 
product of the first equation of (|13h with each of the unit vectors leads to the equations 

4i-n 2 = ij4(t) (16) 
m 

n 1 -n 3 = -i^ 1 (t) . (17) 
m 

Now taking the scalar product of the second equation of (|13j) with each unit vector in turn, and 
using (fl~5l) to (fT7|) to simplify gives, respectively 

^ = ^ + l § 9[R' l2 (t) + R 21 (t)] (18) 

Xf) 1 Xfl 

^ = -[R'22(t) - Rn(t)} - -[R'At) + R' 2l (t)} (19) 

dt> m m 

n 2 n 3 = -R'szit) . (20) 
Using the final equation of f)13[) . and making use of (|15p to (|20p to simplify, we find 

Y^=Y[ + ^084>(z[R' 12 (t) + R' 2l (t)) + # 13 (t) + i4(t)) (21) 

g = l^)-^)]^^^ (22) 

o<p m m m \ J 

Z = + ^3 2 (*)] + " • (23) 

Equations (|15p to (|23p are the exact equations of motion for the nine variables parametrising 
5r. Retaining only the leading-order terms in 50 and 6(f), we have 

~sF~ = ~ R '^ ' Ja = ~^22 ~ R'n) > = IT (#33 ~ R'22) ■ ( 24 ) 

X\ m dv m o<p m 

Using equation ()10p . we find the following expressions for the Lyapunov exponents: 

A: = -(R' u ) , A 2 = -(R'22) , A 3 = -(R 33 ) • (25) 
m m m 

The Lyapunov exponents are thus obtained directly from the elements of (R') which satisfies 
an equation similar to f)12[> : to derive it, write = Oe^ where are the unit vectors of a 
fixed Cartesian coordinate system, were introduced in ([9]), and O is an orthogonal matrix. 
Transforming (|12l) 

O+RO = - 7 R' - — R' 2 + F' . (26) 
m 

Making use of 

O+RO = R' - [R', O+O] (27) 
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we obtain 



R' 



7 R' - — R /2 + [R, O+O] + F' . 



m 



(28) 



The matrix elements of + are given by (|16[) . (|17p . and ([20 



O+O 




R 21 



'32 



#31 \ 



32 





(29) 



To summarise: the Lyapunov exponents are determined by (j25H . using expectation values of 
variables occurring in the system of stochastic differential equations defined by (|28p and (|29p . 



3.2 Numerical calculation of Lyapunov exponents 

Our numerical results for the Lyapunov exponents, represented as symbols in figures [HElEl and 
[6J were obtained in the limit of r — > by discretising time in equation (T7|). We write t = n5t with 
n = 0, 1, 2, . . ., and with a time step 5t which satisfies 1/7 S> <5i 3> r. Writing 5r n = Sr(nSt) 
and <5p n = 6p(nSt) the Euler discretisation of ([7]) is 

£ r ("+i) \ _ / I I(5t/m) 
S p (n+1) J ~ y F^St I(l-7<5t) 

Here F^ n ' is a matrix with elements 

with n, v = 1,2,3, and I is the 3x3 unit matrix. The random components F^ v average to 
zero, are uncorrelated for n / n', and the correlation (Fj£y F^fy) is determined by the statistical 
properties of u(r, t) described in section 14.11 The Lyapunov exponents are identified as the 
asymptotic growth rates of the three largest eigenvalues of the product M^M^ -1 ) • • • M^ ) 
nearly diagonal matrices for large values of N. The asymptotic growth rates are determined 
using a method described in [29J . 



6rH I = mW f 5r[n) 

Sp {n) J — M I Sp (n) 



(30) 



4 Langevin equations for Lyapunov exponents 
4.1 Equations in Langevin form 

The equation (|28p for R' can be simplified when the correlation time of the velocity field is 
sufficiently short, and when the amplitude of the random force is sufficiently small. In this limit 
the force-gradient term F' behaves like a white-noise signal, and the equations of motion reduce 
to a system of Langevin equations: 

dR' = (-7R' - — R' 2 + [R, O+O] V + dC (32) 
V m J 

where d£ is a matrix of random increments satisfying 

(dC^) = , (dC^dC^v) = 2V^, u ,dt . (33) 

The elements 2?„ y „v of the 'diffusion matrix' depend on the statistics of the fluid velocity field 
u(r,t). The latter may be decomposed into potential and solenoidal components, generated 



11 



from scalar and vector potentials. The Stokes drag force on the particle, / = mju is written in 
terms of potentials <\> and A: 

f = + V A A . (34) 

The fields 4>(r, t) = Aq(v, t) and A p (r, t), \i = 1, 2, 3 are assumed to possess statistical properties 
which are homogeneous in space and time, and isotropic in space. Also, it is assumed that these 
fields are uncorrelated, and the intensity of the A p fields (for fi = 1,2,3) is such that the 
correlation function is of the form 

(A^ro + R, t + t)A v (r , to)) = 5 pu [(1 - a 2 )5^ + a 2 ] C{R, t) (35) 

for some constant a, where R = \R\. 

Using spatial and temporal homogeneity of the velocity field, and noting that in the limit 
Ku — > the particle does not move significantly in time r, elements of the diffusion matrix are 

W-i/>(^(0.«>§£<0,0)). (36) 

Inserting the expression (|34p for the force gives an expression for these elements in terms of 
second derivatives of the fields A^. For any isotropic field A(r,t), these satisfy, for example, 

ld 2 A t x d 2 A f \ I d 2 A , . d 2 A . \ Id 2 A, d 2 A \ 

Now consider the evaluation of the elements of the diffusion matrix, Pm/uV- First we note 
that due to isotropy, the values of the elements are invariant under any permutation of indices 
(for example, V 33 ^ 22 = f 11,33) • It is easy to see that V^^i^ = unless the four indices can 
be grouped into two equal pairs. There are four cases where indices are paired, namely T> a a,bbi 
T^ab,ab, T^ab,ba and V aa ^ aa where a and b are different numbers from the set {1, 2, 3}. We define 

(38) 



s„ = i£>(|f(<M)|f(o,o) 



and use (|37p to express the non-zero elements T> fl v,LL'v' m terms of T>o and a 2 . It is simplest 
to calculate specific examples of the non-zero elements, and to deduce the others by permuting 
indices: writing d 2 A p /dr u dr p (r,t) in shorthand as (d 2 p A p ) t , we have 

/oo 
dt ((d 2 u 4> + d 2 21 A 3 - d^AzUd 2 ^ + d 2 21 A 3 - d 2 3l A 2 ) ) 
-00 

f 2a 2 \ 

/oo 
dt ((d 2 2 <f> + d 2 2 A 3 - d 2 2 A 2 ) t (d 2 2 ^ + d 2 22 A 3 - d 2 32 A 2 ) ) 
-00 

(\ 4a 2 \ ^ 

/oo 
dt ((d 2 n 4> + d 2 21 A 3 - d 2 31 A 2 ) t (d 2 2 <p + d\ 2 A x - d\ x A 3 )o) 
-co 



-00 
00 



/oo 
dt ((d 2 2 4> + d 2 2 A 3 - d 2 2 A 2 ) t {d 2 2 4> + d^Ax - d 2 n A 3 ) ) 
-00 

" n '(i-T)- n " (39) 
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where the first three relations define the constants T>\, T>2 and T>%. 

The Langevin equations (|32p can be labelled by a single index, which will be indicated by 
using a Roman letter, and packing the double indices such that i = 3(/i — 1) + v. The above 
Langevin equations are then of the form 



dR' 



1 



9 9 



-r^d* --EE y^R' k 



j=l fc=l 



dt + dd 



(40) 



with (dCid^j) = 2T>ijdt. The Vijk are determined by (|32p ; most of them are zero. It is convenient 
to scale the equations for R', to dimensionless form. We write 



t' = jt 



,1 I 



1 ~zr~ Rj i 



dwi 



Pi s 



and 



V 



1/2 



m7 



3/2 



1 + 4a 2 
3 + 2a 2 ' 



(41) 



(42) 



The parameter e is a dimensionless measure of the particle inertia, and the parameter T char- 
acterises the nature of the flow: it ranges between | and 2, and we have 



for potential flow 



2 for solenoidal flow . 



(43) 



In [19] where two-dimensional flows are discussed, the corresponding parameter satisfies 1/3 < 
r < 3, and solenoidal flow is obtained for T = 3. 
In scaled coordinates, the Langevin equations are 



dxi 



9 9 

+ «EE V ijkXjX k 

j=l k=l 



dt' + dwi = Vidt' + dwi 



where a 'velocity' v with components 



(o) , (i) 



,(°) 



9 9 

■ ^ ^ VijkXjXk 
j=i h=l 



(44) 

(45) 
(46) 



was introduced, and (dwidwj) 
Dij is given by 



2Dijdt'. The transformed diffusion matrix D with elements 





fl 











a 











a \ 









r 





a 



























r 
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a 





r 



















D = 


a 











1 











a 


a = 



















r 





a 















a 











r 



























a 





r 









U 











a 











iy 




Note that T = V 2 /Vi and 


a 


= v 3 /v 


L- 


The diffusion matrix 



Ci-r) 



(47) 



form, with three 2x2 blocks and one 3x3 block, which is in cyclic form: it can therefore be 
diagonalised analytically. 
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4.2 Discussion 



The problem of determining the Lyapunov exponent of particles suspended in a turbulent fluid 
has thus been transformed into the problem of determining expectation values of the position of 
a particle undergoing a Langevin process in a nine-dimensional space (defined by equation il 1 II:;. 
From (|4ip we see that {R'^/m = je(xi), and from (|25p we see that the Lyapunov exponents are 

Ai = 7e(zi> , A 2 = je(x 5 ) , A 3 = -/e(x 9 ) . (48) 

The velocity ()45|) in the Langevin equation contains terms which are linear in the displacement, 
D", which drive this Langevin particle back towards the origin. If, however, the particle diffuses 
sufficiently far from the origin, the quadratic terms t/ 1 ) in the velocity become dominant, and 
may take the particle away to infinity. In fact, numerical studies show that the particle always 
does escape to infinity. We must consider the significance of this effect. When the inertia of 
the suspended particles is large, the momenta of the particles are not determined solely by their 
positions. It is therefore possible for two of the vectors 5r^ to become coplanar, while the 
vectors 5p^ continue to span three dimensions. As the vectors 5r^ approach co-planarity, the 
inverse of the matrix R defined by (llip must become singular. Correspondingly, some or all 
of the components Xi must diverge to infinity in a finite time. After the point x{t) diverges to 
infinity, continuity of 5r^ and 5p^ implies that it immediately reappears and converges from 
the reflected point at infinity. In terms of the parametrisation of the displacements 5r fl given 
in ([9]), we see that two of the vectors become collinear when 5(j) = and all three are collinear 
when 56 = 0. The 5(f> = is therefore a co-dimension one-condition, which is realised by varying 
time. When 6(f) = the elements of the second and third columns of R' diverge. Generically, 
the condition 59 = never occurs. 

The events where the elements of R' diverge in the Langevin simulation therefore correspond 
to events where the local density of the suspended particles diverges because the volume element 
of their flow vanishes [20]. These events are termed caustics. Simulation of the Langevin 
equations therefore also gives the rate J at which an particle passes through a caustic (which is 
identical to the rate at which Langevin trajectories escape to infinity), as well as an estimate of 
the Lyapunov exponents. 

The Langevin equations (|32p are a valid approximation for (|28p and (|29p provided that two 
conditions are satisfied. The correlation time of the forcing terms must be small compared to 
the relaxation time 7 : this clearly implies that the Stokes number should be large, that is 
0<1. A second condition is that the random forcing term should be sufficiently weak that the 
displacement of the coefficients R^ v during the correlation time r should be small, relative to 
their typical values. This condition was discussed in [21J: in the notation of this present paper 
it leads to the requirement that k -C 1. The arguments of this section are therefore valid in the 
limit where SJ< 1 and n -C 1. 

Figure [3] illustrates the application of the Langevin equations (|44j) to calculating the Lya- 
punov exponents for a flow with a very short correlation time. The results are compared with 
a direct evaluation obtained by multiplying the monodromy matrices of the corresponding flow, 
as explained in section 3.2. The first three panels show the Lyapunov exponents as a function of 
e for three different values of T, corresponding to incompressible flow (r = 2), purely potential 
flow (r = |), and a mixed case, T = 1. The final panel illustrates the limiting behaviour of the 
Lyapunov exponents as e —* 00. In that limit the effect of the damping 7 becomes negligible, 
and the Lyapunov exponents reach limiting values that are independent of 7. This implies that 
the functions fj(e) = Aj/7 are asymptotic to 

m = k „ (49) 
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in the limit as e — ► oo, for some coefficients Cj which depend upon T. We have not been able to 
determine the coefficients Cj(T) analytically. 

The Langevin equation does not appear to be exactly solvable. Section [5] considers how to 
obtain a perturbation series solution for the variables (a;,). 

4.3 Estimate for diffusion constant in turbulent flow 

In section 14.11 we discussed the definition of e for a flow with small Kubo number. However, 
the principal area of application of our results is to particles suspended in a fully developed 
turbulent flow. In section [2~?2l we argued that the Kubo number is always of order unity for this 
case. Here we consider how the definition of e must be adapted to make our results applicable 
to turbulent flows, and what is the appropriate value for the correlation time r. 

The dimensionless parameter e is expressed (equation (|42p ) in terms of a diffusion constant 
defined by (|38|) . This quantity is related to the spectral intensity of the rate of strain in the 
neighbourhood of a particle with trajectory r(t), defined by 

I(u) = dt exp(i^)^(r(t),i)^i(O,0)^ . (50) 

Note that (unlike (|38p ) we consider the temporal variation of the position r(t) of the suspended 
particle, because when Ku = 0(1) this may change by a significant amount during the correlation 
time, r. Let us first consider how to estimate this quantity when the Stokes number of the 
suspended particles is small, so that the correlation function in (|50p may be approximated by a 
Lagrangian correlation function. We estimate T> oc m 2 7 2 I(0), so that e 2 oc 1(0) /j. The intensity 
1(0) has the dimension of inverse time, so that 1(0) oc t~ , where t s is a characteristic time 
scale of the flow. In the case of fully developed turbulence, it is not immediately clear whether 
r s is the time scale associated with the dissipation scale, or else some longer time scale. The 
Kolmogorov theory of turbulence [26j implies that if the integral defining I(oS) is dominated by 
the inertial range, we may write /(cj) in terms of £ , the rate of dissipation per unit mass, with 
no dependence upon v. We therefore seek a relation of the form 

/(cj) = CrV 3 (51) 

for some constants a, (3, C. This relation is only dimensionally consistent for a = 0, f3 = 1, 
which makes the result independent upon £ (and furthermore implies that it vanishes at cj = 0). 
We infer that the value of T>q is not determined by the inertial range of turbulent flow, and we 
expect that it is determined by the dissipative scale, with characteristic time scale r ~ \J u/£. 
We therefore expect that 

(52) 

7 V v 

for fully-developed turbulence, where Kolmogorov's 1941 theory of turbulence [26J predicts that 
K is a universal constant. 

Finally we consider two refinements of the estimate (152p . Recent insights concerning inter- 
mittency suggest [26J that K should have a weak dependence upon Reynolds number of the 
turbulence. More significantly, at very large values of the Stokes number, equation (|52p may 
overestimate e 2 since d Xl ui(r(t), t) will fluctuate more rapidly than it would for a particle which 
is advected. (An analogous effect is described in detail in [30].) These considerations suggest 
that we write 

2 _ if (St, Re) 

7 
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where St = \JE jv^ 1 is the Stokes number referred to the Kolmogorov timescale. The function 
-fT(St, Re) approaches a finite limit as St — ► 0, but approaches zero as St — ► oo, and has a weak 
dependence upon the Reynolds number, Re. 



5 Perturbation theory 

5.1 Mapping to Hamiltonian form 

The Langevin equations (|44j) are equivalent to a Fokker-Planck equation for the probability 
density P(x,t f ) of x 

dP 

_ = V-(-u + DV)P (54) 

When e = 0, the velocity v of the Fokker-Planck equation (|54p is linear in the displacement from 
the origin, see equations (|45l46p . It is easily verified that the solution is a Gaussian function. 
This suggests that it may be possible to map the unperturbed (e = 0) problem to a nine- 
dimensional harmonic oscillator. We shall show how this mapping can be achieved, and what 
is the form of the perturbation representing the velocity terms which are quadratic functions of 
the coordinates. 

We write the Fokker-Planck equation as 

dP 

— = FP = (F + eF 1 )P (55) 

where the notation A indicates that the object A is an operator. We have 

F = V • (-v {0) + DV) , Fi = -V ■ (56) 



with 

„(0) = _ x 

and where the components of t)W are the quadratic terms in v proportional to 
e, defined in (|46p . We consider the steady-state solution which solves FP(x) = for e = 0: this 
is 

P (x) = ^exp(-ia: • D 1 x) = Aexp[-$ (x)] (57) 
(where A is a normalisation constant). The latter identity defines &o(x). We now define H by 

H = exp($o/2)Fexp(-$ /2) (58) 

and re-write the steady-state Fokker-Planck equation as 

HQ{x) = (59) 

where Q(x) = exp[<&o(x)/2]P(x). This form will be suitable for a perturbation expansion in 
the small parameter e, since H consists of two parts, i.e., H = Hq + eH±, with Hq a Hermitian 
operator. 

We require an explicit expression for Hq = exp($o/2)-^b ex P( — ^o/2)- By considering the 
action of Hq upon an arbitrary function / we find 

= U2^u ~ l^2x i (D' 1 ) ij x j D ij d?j 

i ij ij 

= itr(I) - \x ■ D- 1 x + d x ■ D d x (60) 

where I is the 9x9 unit matrix and d x = (d Xl , . . . , d Xg ) = (d\, . . . , dg). Note that ([UUP is the 
Hamiltonian operator of an isotropic nine-dimensional quantum harmonic oscillator. 
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Next consider Hi = exp($>o/2)Fi exp(— &q/2): again we consider the action of H\ upon an 
arbitrary function / 

-Htf = ^exp(cD /2)5 i K (1) exp(-$ /2)/) 

i 

= J2 d M 1] f) -ffirt 1 // (si) 

i ij 

so that 

H! = -d x ■ »W + • D ^cc . (62) 
5.2 Use of annihilation and creation operators 

We next consider how to re-write the operators Hq and H\ in terms of harmonic-oscillator 
creation an annihilation operators. This is easier if we first diagonalise Hq. The diffusion matrix 
D is a real symmetric matrix, which can be diagonalised by an orthogonal matrix, U 

U _1 DU = A (63) 

where A is a diagonal matrix with elements Ay = Aj<5y, and Aj are the eigenvalues of D. We 
find 

H = ±tr(I) - \ V ■ A" 1 rj + d v ■ A d v (64) 

where 

T) = U 1 x , d v = IT 1 d x . (65) 
Since A is diagonal, equation (j64"j) simplifies to 



^o = l + |E(2A^-^f). (66) 



Now define the following operators 



/ d 

Pi = -iV2Ai — 
drji 

* = ^m* (67) 

which are analogous to the position and momentum operators of quantum theory. In terms of 
qi and pi (|66|) becomes 

i 

where / is the identity operator. This is the Hamiltonian of an isotropic quantum harmonic 
oscillator in nine dimensions. The commutator of the qi and pi operators is [<ji,Pj] = HSij, so 
that the pairs qi and pi have the same commutators as the position and momentum operators 
in quantum mechanics (with h = 1). We next introduce 

Qi + i Pi 



V2 
qi-ipi 



<69) 
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whose commutator is [dj, ctj ) = The af and dj are respectively the creation and annihilation 
operators, or raising and lowering operators, for the degree of freedom labelled by i. Now we 
find 

9 

H = -J2afai- ( 7 °) 

i=i 

Reducing the unperturbed operator to this form is useful because the simple algebraic proper- 
ties of the annihilation and creation operators make it possible to perform perturbation theory 
exactly, to any desired order. This is achieved by using the set of eigenfunctions of Hq as a conve- 
nient basis set. These eigenfunctions are labelled by a set of quantum numbers {mi,m2, ..,7719}, 
taking values rrij = 0,1,2,3.... Thus each eigenfunction is labelled by a vector m with non- 
negative integer elements. The eigenfunction with this set of quantum numbers will be denoted 
by a Dirac 'kef vector, \m) [31J. (We use this notation rather than the more usual \m) to 
avoid confusion with the angular brackets denoting averages.) For the Hamiltonian (|70p . the 
eigenvalues are just the sum of the quantum numbers, so that the eigenvalue equation is written 

9 

Ho\m) = — >Jmj|m) . (71) 

The annihilation dj and creation af operators map one eigenstate to another, by (respectively) 
raising and lowering the quantum number mf 



aj\mi,rri2, --jTrij, ..,mg) = J rrij \mi, ni2, .., rrij — 1, .., mg) 



af \mi,m 2 , ..,rrij, ..,m 9 ) = Jmj + 1 |mi,m 2 , ..,m,j + 1, ..,m 9 ) (72) 

and if a,j acts on an eigenstate for which the quantum number rrij is zero, the result is zero. 

We now consider how to rewrite H\ in terms of the dj and af operators following the same 
procedure as for Hq. Using ([65]) we have 



d x ■ = J2 UAvf , D 1 x = ]T v^Uji Kfm . (73) 



,f , r -I) .r X 

10 ij 

We obtain 

Hi = -E^ " kK\) «i W (*(»»)) (74) 
ij v h 7 

where we have made the dependence of the components of v^ 1 ' (x(r])) upon the variables r] 
explicit. 

The final step is to express H\ in terms of the dj and af operators. Inserting rji = (ai+af)\fKi 
and d/drji = (dj — af)/2\/X i into equation ([71]) yields 

H 1 = J2U ll ^=afvf(x( V )) . (75) 

ij * * 

Now each of the -uj 1 ^ can be expressed as a quadratic form 

vf> (x) = - E aw (76) 

where Vj/y are the coefficients (almost all zero) of the terms x^xi appearing in each . Since 
03(77) = U77 and rji are given by rearranging (|69|) . we have 

x fc = E u km{a m + d+)v / A^' • (77) 
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Combining with (I76p and inserting into (I75p gives 

Hi = £ #2* «* + («™ + + &+) (78) 

imn 

where 




Note that the coefficients H^ n defining the perturbation operator Hi can be obtained exactly, 
because the matrix D can be diagonalised exactly. 



6 Iterative development of the perturbation series 
6.1 Method for developing the series expansion 

Instead of solving the Fokker-Planck equation F P(x) = we attempt to solve H Q(x) = 0, 
where Q(x) = exp[&o(x)/2]P(x). In the following we use a shorthand notation for integrals 
which is equivalent to the Dirac notation [31J of quantum mechanics: given two functions a{x) 
and b(x) and an operator A, we define 

/poo poo 
dx a*(x) Ab(x) = / dx\ . . . dxg a*(xi, . . . , xg) A b(x±, . . . , xg) . (80) 
J — CO J — oo 

Here the asterisk denotes complex conjugation. Now consider how to obtain the Lyapunov 
exponents from the function Q. They are obtained from the mean values (xj) of the coordinates. 
These are 

(x^ = JdxxiP(x) 

Jdx exp[-<& (x)/2}xiQ(x) 

J dx exp[-$ (x)]Q(x) 
(Ojgijg) 

(o\Q) 

= j^J2 u ^(°\^ + ^\Q) (81) 

where we have used the fact that exp[— &q(x)/2] is the eigenfunction of the 'ground state', |0). 
The denominator is included to take account of normalisation. 
We calculate \Q) by perturbation theory: writing 

\Q) = \Qo)+e\Q 1 )+e 2 \Q 2 ) + ... (82) 

we find that the functions \Qk) satisfy the recursion relation 

\Q k+1 ) = -H l Hi\Q k ) . (83) 

At first sight this appears to be ill-defined because one of the eigenvalues of Hq is zero, so 
that the inverse of Hq is only defined for the subspace of states which are orthogonal to the 
'ground state', |0). However, because all of the components of H\ have a creation operator as 
a left factor, the state is orthogonal to |0) for any function \ip), so that ([83]) is in fact 

well-defined. The iteration starts with \Qq) = |0). 
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The functions \Qk) are expanded in terms of harmonic-oscillator eigenstates |m), with coef- 
ficients dm and m = (mi, . . . , mg): 

\Q k ) = J2^\m). (84) 

m 

By projecting equation (|84p onto the vector |m) and using the fact that the eigenvectors \m') 
of Hq form a complete basis, the iteration can be expressed as follows (for m ^ 0): 



„(fc+l) _ V i m \ H l\ m ') Ak) /qk\ 



The matrix elements {m\H\\m') are readily computed using the algebraic properties of the 
raising and lowering operators, ([72]) [31]. The coefficients a™ are then calculated recursively, so 
that we obtain the states \Qk)- Finally, expectation values are extracted using (I81j) . 

6.2 Programming the perturbation expansion 

It is not practicable to carry out the perturbation expansion by hand even at leading order, 
because the Hamiltonian H\ contains several hundred non-zero coefficients H^jj,. The calculation 
was automated using a Mathematica program. 

From (I83h we see that the n th order of the perturbation expansion requires the calculation 
of 'matrix elements' of the form 

I(n,j) = (OKoj- + a+)F - 1 i? 1 F - 1 Fi . . . H^H^O) (86) 

with Hi occurring n times. Because the Hamiltonian has been expressed in terms of raising and 
lowering operators, each successive application of Hi or Hq 1 to the eigenfunction |0) gives a 
function which consists of a linear combination of a finite number of eigenfunctions |m). The 
inner product is then readily evaluated using the orthonormality of the eigenfunctions: 

9 

(m\m f ) = b m , m , = Y[ £ mi ,m' . (87) 



The number of eigenfunctions included in the sum (|84p increases very rapidly as the order of 
the perturbation increases. The matrix element I(n,j) can be written as the inner product of 
two function vectors: 

= (x(k,j)\ip(n- k)) 
\ X (k,j)) = HiH 1 HiH 1 ...HiH 1 (a 1 +a+)\0) 

v 

.ffi appears k times 

|^(m)) = H a l Hi...H^ l Hi\Q) (88) 

" v ' 

Hi appears m times 

where the choice of k is, in principle, arbitrary. However, the computational effort is proportional 
to the sum of the number of terms in the vectors \x( n , k)) and \ip(m)), and this is minimised by 
taking fc k n/2 in equation ([88]) . 

The program was checked by comparing the coefficients for Ai with those determined in [21] 
using a perturbation of a two-dimensional harmonic oscillator, which allowed evaluation of the 
first 47 non-vanishing coefficients. We also remark that for the second Lyapunov exponent it is 
sufficient to consider a system of six coupled harmonic oscillators, because the Langevin equa- 
tions for (xi, . . . , xq) do not depend upon the values of (X7, x$, xg). This allows the perturbation 
series to be taken to higher order for the second Lyapunov exponent. 
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6.3 Results 

We find that the denominator (0|Q) in (|8ip is equal to unity to all orders in e, and from the 
numerator we obtain series expansions for the Lyapunov exponents in the form 

oo 

V7 = E^OV"- (89) 

n=l 

Note that only even orders in e occur in this expansion: this fact is most readily understood 
using an argument which will be given in section [71 The first seven coefficients Cn are 

c^\r) = -1 + 2T 

c W(r) = -5 + 2or-i6r 2 

4 1} (r) = -60 + 360r- 568r 2 + 272r 3 

c W(r) = -1105 + 8840 T- 61936 T 2 /3 + 58432 T 3 /3 - 19648 T 4 /3 
C W (r) = -27120 + 271200 T - 7507040 r 2 /9 + 3492160 T 3 /3 

-2316032 T 4 /3 + 1785856 T 5 /9 
4 1J (r) = -828250 + 9939000 T- 1020068800 T 2 /27 + 1874157440 T 3 /27 

-613664384 T 4 /9 + 934756352 T 5 /27 - 193558528 T 6 /27 
4 1} (r) = -30220800 + 423091200 T - 154727293760 T 2 /81 + 351319669760 T 3 /81 

-454943581760 T 4 /81 + 342675611776 T 5 /81 - 140392616704 6 /81 

+24271797760 T 7 /81 

4 2) (r) = -2 + r 

4 2) (r) = -i6 + 22r-5r 2 

4 2) (r) = -549/2 + 1287T/2 - 871 T 2 /2 + 125 T 3 /2 

4 2) (r) = -13463/2 + 22506 T- 79121 T 2 /3 + 35606 T 3 /3 - 7723 T 4 /6 (90) 
4 2) (r) = -627719/3 + 2731795 T/3 - 13765330 T 2 /9 + 3598630 T 3 /3 

-1229363 T 4 /3 + 341381 T 5 /9 
4 2) (r) = -280669811/36 + 25085281ir/6 - 9891631295 T 2 /108 + 2783144725 T 3 /27 

-2202644629 T 4 /36 + 924868595 T 5 /54 - 157226321 r 6 /108 

4 2) (r) = -145680639449/432 + 928376776943 T/432 - 7524732877927 T 2 /1296 

+11046167913985 T 3 /1296 - 9388185321985 T 4 /1296 + 4517559789671 T 5 /1296 
-1087174658765 r 6 /1296 + 87859310987r 7 /1296 

c f)(D = -3 

4 3 ) (r) = -33 + 12T 

4 3) (r) = -1479/2 + 1215T/2 - 429 r 2 /2-3T 3 /2 
4 3) (r) = -45627/2 + 29954 T- 21071 T 2 + 5394 T 3 + 143 T 4 /2 
4 3) (r) = -1731931/2 + 18888175 T/12 - 4963475 T 2 /3 + 5130835 T 3 /6 
-1013585 T 4 /6 - 37765 T 5 /12 
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c [j 3) (r) = -461447213/12 + 268463707T/3 - 2249372585 T 2 /18 + 872716325 T 3 /9 
-1414393915 T 4 /36 + 56257766 T 5 /9 + 1316179 T 6 /9 

4 3) (r) = -279693223943/144 + 1190601381865 T/216 -4142788537873, T 2 /432 

+2141088699035 T 3 /216 - 290883665975 T 4 /48 + 144871714325 T 5 /72 
-38215846457 T 6 /144 - 1592440259 T 7 /216 . 



Equations (|89p and (190p are the main results of this paper. They provide an expansion of the 
Lyapunov exponents in terms of the dimensionless parameter e ~ k£1~ 1 / 2 . The expansion is valid 
when k <^ 1 (small Kubo number) and 0^1 (large Stokes number; this is the underdamped 
limit where inertial effects are dominant). 

The coefficients of V in (j90|) are all rational numbers, despite the fact that the algebra of the 
raising and lowering operators (equations (|72|) ) produces expressions containing square roots of 
integers. The reason for cancellation of the square roots to produce rational coefficients is not 
yet understood. 

In the remaining sections we discuss physical and mathematical implications of (|89p and 
(|90l) . and relate these results to results obtained in the overdamped limit (where inertial effects 
are small). 

6.4 Summation of the perturbation series 

The coefficients in (|90p grow rapidly with order n, indicating that these are divergent series. We 
find they have factorial growth for large n: 

c^~C5 n (n-l)! (91) 

for some constant C. This factorial growth is typical of an asymptotic series [33]. We find that 
S = 1/6, independent of j and T. Figure H] shows the growth of ch^ for T = 2. 

There are rather few concrete physical problems for which perturbation expansion coefficients 
are available to high order: most studies are concerned with quantum-mechanical problems which 
are perturbations of a harmonic oscillator, such as [32]. Here we do have the perturbation series 
coefficients to high order, for the same underlying reason: we used the algebra of harmonic- 
oscillator raising and lowering operators to compute matrix elements exactly. 

Methods for treating divergent series are discussed in [33] and [34] , assuming that the expan- 
sion can be pursued to high order. Here we use Pade-Borel summation, similar to an approach 
used to sum the perturbation series of certain one-dimensional quantum-mechanical anharmonic 
oscillators |32j . We evaluate the 'Borel sum' 

^ma x (J) 

= E V" (92) 

i lb. 

71=1 

where n max is the number of terms available. The sum of the series is then taken to be 

A i / 7 = Re J dtBj^t)^ 1 (93) 

where C is a suitable curve in the complex plane. Assuming that the Borel sum has a finite 
radius of convergence, the function may be approximated by Pade approximants of order N/M, 
PN/M{ e2 t) [35]. Here iV and M are the orders of the numerator and denominator polynomials, 
respectively and N + M < n max , where n max is the number of coefficients available. 
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The integral in (j93[) is evaluated numerically. If the Pade approximants to Bj(e 2 t) have poles 
on the positive real axis, the integration path in (|93p is taken to be a ray in the upper right 
quadrant in the complex plane. 

Results of Pade-Borel summation of the series for Ai in the incompressible case (T = 2) 
are shown in figure [5k for N = M = 5, 10, 15, 20, 23. The summation results are in good 
agreement with numerical results for Ai obtained as described in section 13.21 For the first 
Lyapunov exponent, we obtained higher-order coefficients using the method described in [21] 
(which is particularly efficient but also restricted to calculating the maximal exponent): in this 
case n max = 47. We conclude that the Pade-Borel approach works very well for the largest 
Lyapunov exponent. 

Figure 03 shows results for the first and second Lyapunov exponents for T = 1, representing 
a flow field with both solenoidal and compressible components. For the maximal Lyapunov 
exponent, Ai, Pade-Borel summation gives very good agreement with the numerical results. 
For the second Lyapunov exponent, inspection of the coefficients (|90l) shows that for T = 1, 
C; = — cl i , in other words, the perturbation coefficients for Ai + A2 vanish identically for 
r = 1. However direct numerical simulations show that Ai + A2 is not equal to zero when T = 1. 
The WKB analysis summarised in section [7] gives rise to the hypothesis that there may be a 
non-analytical contribution to the Lyapunov exponents of the form Cj(F, e) exp(— l/6e 2 ) which 
has to be added to the Pade-Borel approximation. Figure 03 shows that this is indeed the case 
for the second Lyapunov exponent for T = 1: adding C exp(— l/6e 2 ) (with C = 0.79 a fitted 
constant) to the Pade-Borel sum for A2 (solid line) gives very good agreement with the numerical 
data. This observation shows that there are contributions to the Lyapunov exponents which are 
not captured by a simple application of the Borel-Pade summation technique. 

Figure [6] shows a comparison of the results of Borel-Pade summation with Monte-Carlo 
simulations of all three of the three Lyapunov exponents for T = 2 (incompressible flow) . Here 
we used the coefficients given in section 16.21 with n max = 7. The results for the first Lyapunov 
exponent show excellent agreement, as mentioned above. Those for the second exponent show 
a systematic deviation for larger values of e. The results for the third exponent show some 
instability upon changing the number of terms included in the Borel sum, implying that more 
terms are required to achieve a convergent result. 

We conclude that the Borel-Pade summation gives very satisfactory results for small values 
of e, but for large e there are systematic deviations which are not yet understood. These appear 
to be associated with non-analytic contributions of the form exp(— $*/e 2 ). 



6.5 Relation to clustering 

A major motivation for calculating the Lyapunov exponents is to evaluate the dimension deficit 
A (defined by equation (JT])) and hence the Lyapunov dimension, = 3 — A. To lowest order 
in e, equations (|89|) and (f90l) imply that 

Ai/ 7 = (-l + 2r)e 2 , A 2 / 7 = (-2 + r)e 2 , A 3 / 7 = -3e 2 . (94) 

In figure [3] a, b, and c these expressions are shown as dashed lines. Equations (|94p imply that 
in for an incompressible flow (r = 2), the sum Ai + A2 + A3 vanishes in the limit of e — ► 0, and so 
does A, defined in (pQ). This reflects that in the absence of inertial effects, an incompressible flow 
cannot give rise to density fluctuations. The leading-order behaviour of the dimension deficit is 

A = 3^6 + 54-54r +i 2 ir% 2 + 0(e4) ^ 
In the incompressible case this gives A = 10 e 2 for small e. 
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Figure shows numerical results for A as a function of e for incompressible flow, including 
the asymptotic approximation (|95p . When the numerical results were re-plotted in figure QJ>, 
we used the relation St = e 2 /0.256 2 , with the factor chosen to give a good agreement as judged 
by eye. 

Figure ©J shows results for the dimension deficit A obtained from Pade-Borel summations 
for the Lyapunov exponents shown in figure [6]a, using equation ([I]). For small values of e the 
agreement is excellent. Despite shortcomings of the Borel summations for A2 and A3 mentioned 
above, we observe satisfactory agreement between the theory and results of computer simulations 
for all values of e where A is positive. In summary, Pade-Borel summation of the series (|90p 
provides a satisfactory theoretical description of the dimension deficit A. The maximal value of 
A obtained here (and in [3]) is in good agreement with direct numerical simulations of particles 
suspended in a Navier-Stokes flow (see figure [Q). 



7 WKB analysis 

In this section we discuss the stationary solution of the Fokker-Planck equation (|54[) by means of a 
WKB method, along the lines discussed in [36] . This gives an indication as to the interpretation 
of the non-analytic contributions to the Lyapunov exponents of the form exp(— <3?/e 2 ) which 
were considered in section [6] A precise evaluation of these contributions is beyond the scope of 
existing methods. 

It is convenient to re-scale the variables Xi in (|54h using x\ = x%e. The stationary Fokker- 
Planck equation for P(x') is then 

V • [(-«' + e 2 DV')iV)] = ( 96 ) 

with v\ = —x\ — J2jk Vijkx'jx'k- Equation ([96]) indicates why the odd orders of the perturbation 
series developed in section [6] vanish: because this equation contains e 2 rather than e, expansion 
of any quantity as a power series must be a series in powers of e 2 . 

In the remainder of this section we drop the primes and make the WKB ansatz [35 1 



P(x) = Aexp 



(97) 



(where A is a normalisation factor). Our objective is to determine the nature of the function 
in (|97p . Inserting this into (|96p and keeping only the lowest-order terms in e gives 

v ■ V$ + V$ • DV$ = . (98) 

This has the form of a Hamilton- Jacobi equation [36j. Recall that given a Hamiltonian H(x,p) 
the Hamilton- Jacobi equation for the action &(x,E) is 

H(x,V$) = E. (99) 

The Hamilton-Jacobi equation is solved by finding the classical trajectories which are solutions 
of Hamilton's equations of motion. The function <&(x,E) is then obtained by integration along 
the classical trajectory 

rt(x,E) 



$(x,E) = dtx-p. (100) 

Jo 

In our case, E = and the Hamiltonian is 

H(x,p) = v(x) p + p-Bp (101) 
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and (using the fact that D is a symmetric matrix) Hamilton's equations of motion are 



dH dH 
x = — =v(x)+2Bp, p = -— = -Fp (102) 

where the matrix F has components i*y = dvi/dxj. Close to the origin, $(a;) is expected to be 
approximately &o(x) = ^x T)x, in agreement with the Gaussian approximate solution (I57jl , 
The appropriate initial conditions for the trajectories are therefore infinitesimally close to the 
origin in cc-p-space with initial condition p = V<I>(;e) = x. 

Singular points of &(x) are points x* where V<J> = 0. Consider a trajectory from the origin 
passing through a singular point x*. At this point vanishes and consequently also p. From 
this point onwards, p = p = and after reaching the singular point the dynamics is thus 
advective, x = v(x), with remaining constant at <£* = $>(x*). Singular points are expected 
to give rise to non-analytic behaviour of the Lyapunov exponents of the form exp(— <E>*/e 2 ). 

There is a trajectory of equations (|102j) for which all of the coordinates except x\ are zero, 
for which there is a singular point at x\ = —1 with action <£* = &(x*) = 1/6. We can therefore 
expect non-analytic contributions to the Lyapunov exponents of the form F exp(— l/6e 2 ), where 
F may have an algebraic dependence upon e. 

The formation of caustics is associated with escape of the Langevin trajectory to infinity. The 
Langevin equations have an attractive fixed point at x = 0, but particles that diffuse sufficiently 
far from this fixed point will escape to infinity. We expect that the escape rate contains a factor 
exp[— <3>*/e 2 ], where 3>* is the action of a trajectory to a singular point. We find that the rate of 
caustic formation is of the form 

J'~J^exp(-<^ austic /e 2 ) (103) 

Numerical results confirming this expectation are illustrated in figure which shows that the 
action associated with the formation of caustics is ^caustic ~ 0.125. 

8 Advective approximation 

In the remaining two sections we discuss the overdamped limit 1! > 1 and consider how this 
connects with our results for the underdamped limit described in sections 3 to 7. A surprising 
fact is that the leading-order formulae for the expansion of the Lyapunov exponents in n are 
identical in the limits f2 ^ 1 and fi> 1 (although the higher-order terms differ). When f2 ^ 1 
the particles are advected by the flow. However, we show that when Sl< 1, despite the fact that 
the particles are not advected, the equations determining their separation are identical to those 
of advected particles in the limit e — * 0. It is this advective approximation which is discussed in 
this section. 

In the advective limit it is of interest to compute the sum of the three Lyapunov exponents 
for particles suspended in an incompressible fluid beyond the leading-order approximation. This 
case is important because its sign determines whether particles in an incompressible flow cluster. 
The result for purely advected particles is zero, so that a theory going beyond the advective 
approximation is required. This special case is considered in section [9l 

This section is structured as follows. First in section [HTTl we discuss the advective approxi- 
mation and the conditions for its validity. We also show that the Lyapunov exponents can be 
obtained from an ltd type Brownian process. The Lyapunov exponents of this process were 
calculated by Le Jan [13j . using a notation and terminology which are difficult to relate to our 
own. In section 18.21 we calculate the Lyapunov exponents using a simpler method, expressed 
in our own notation. We calculate the leading-order expressions for the Lyapunov exponents 
for the model with force statistics defined by equation (|35p . considered as an expansion in k, 
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showing that these are identical when £1 3> 1 to those expressions obtained in section [6] when 
0< 1. 



8.1 Brownian advection model 

Consider the linearised equations of motion. These can be integrated to give 

-At 



m Jo 

and 

rti 



1 

Sr(At) - 5r{0) = — / dii 5p(h) (104) 
m Jo 



Spfa) = I 1 dt 2 exp[-7(ti -t 2 )]F(r(t 2 ),t 2 )Sr(t 2 ) . (105) 

J — oo 

We now consider the mean and variance of the displacement (|104p . 
Mean value of displacement 

Combining these expressions, writing the result in component form, and expanding the functions 
in the second integral about r(0), we obtain 

Sr^At) - Sr^O) = — VW dti / dt 2 exp[- 7 (ti - t 2 )) F^(r(0), t 2 )Sr v (0) 

m^Jo J-oo 

+^2J / d M d*2 exp[-7(ti - tjj)] / dt 3 / dt 4 exp[- 7 (t 3 - 1 4 )] 
m z ^ Jo J-oo Jo J-oo 



92 U -(r(0),t 2 )f x (r(0),t 4 ) + §^(r(0),i 2 ) |^(r(0),t 4 ) 



5r v (0)5r A (0) 



ar v ar\ or v dry 

+0(f 3 ) . (106) 

This approximation is valid when the absolute displacement Ar during time At is small com- 
pared to the correlation length: |Ar(At)| <C The force is derived from potentials, as specified 
by equation (|34p . The potentials <ft and have statistics (given by (|35p ) which are spatially 
homogeneous and isotropic. Consider the correlation function of two such quantities A and B 
which have spatially homogeneous statistics. The following relation holds between correlation 
functions involving derivatives: 

^(r,t)B(r',t)+A(r,t)^(r',t'^ = . (107) 

Applying this relation to (|106p . with df fJi /dr u and f u playing the role of the functions A and B, 
we see that 

(Sr(At) - <5r(0)) = 0{At 2 ) . (108) 

Variance of displacement 

The variance of (|104p is 

1 f°° 

(\Sr(At) - 5r(0)] 2 ) = At— / dt (Spit) ■ Sp(0)) + OiAt 2 ) . (109) 

m z J-oo 

The correlation function of the components of the momentum difference is 
<%(<)<$jy(0)> ~ dt i / di 2 exp[ 7 (ti+t 2 )] 

vv , J-oo J-oo 

x ^(rCt + tiJ.t + tiJ^CrCtaJ.tajNtfr^O^vCO) . (110) 
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We consider the evaluation of the momentum correlation function (IllOp in successively the 
underdamped and overdamped limits. In both cases we assume that the displacement Ar 
during the relevant correlation time is small compared to £. In the underdamped case, f2 ^ 1, 
we can approximate the correlation function appearing in (jllOp by a delta function multiplied 
by a weight 

^(r(ti),tx)^(r(t 2 ),t 2 )\ ~22V )MV S(ii-i 2 ) (111) 

where the coefficients V^y^ are defined by equation ([33]) . Using this approximation to simplify 
(fTTUj) we find 

(<5p M (t)<5jy (0)) = -exp^liD^P^^v £rv(0)5rv(0) , (underdamped case), < 1 . (112) 

In the case where 7T » 1, the integral (jllOp is dominated by the region t\ ~ t 2 ~ 0, and we 
obtain 

(<5p M (*)%'(())) = J] ^(^^^'^^^(O)' )) M0)*»v(0) , (overdamped case), fi » 1 . 

(H3) 

Thus we find that the momentum correlation functions are different in the underdamped and 
overdamped cases. However, when we evaluate the variance of the change in spatial separation 
using (|109p . we find the same result in both limits: 

2At 

([Sr(At) - 5r(0)] 2 ) = ]T V„, 5r u (0)5r u ,(0) . (114) 

tyi 7 . 

Equations (I108p and (I114p give the statistics of the change in the particle separation after a 
short time At. This result is valid both when $7^1 (overdamped limit) and when !] C 1 
(underdamped limit), provided that the particle displacement is sufficiently small . In both 
cases we estimate the displacement as |Ar(At)| ~ V DAt, where D ~ u 2 t is the spatial diffusion 
constant. In the underdamped case the correlation time is 0(1/7), and the requirement that 
Ar(7 _1 ) £ gives the condition u^/r/^^/7 1 that is e <C 1. In the overdamped limit (|108|) 
and (I114p are always valid when k <C 1. 



Relation to advective model 

Equations (|108p and (|114p imply that both in the overdamped limit and in the underdamped 
limit when e < 1, the dynamics may be approximated by a random advection model. In this 
model the displacement of a particle during time At is a function only of the position of the 
particle at time t, and not its momentum: we write the displacement at t = nAt as 

r(t + At) = r(t) +w n {r)At (115) 

The values of w n at successive time steps are uncorrelated and average to zero, but w n (r) 
is a smoothly varying function of position, so that nearby particles follow similar paths in this 
random advection model. The spatial correlation of w n (r) at equal times (n = n') is determined 
by that of u(r, t). 

The separation 5r of two trajectories satisfies 

5r(t + At) = [I + M n ]Sr(t) (116) 



27 



where I is the identity matrix and where M n is a matrix with elements (M,,)^ = {d{w n ) il / dr v )5t. 
The elements of this matrix are assumed to have mean value zero and to have diffusive incre- 
ments: for consistency with (I114p we write 

((M n )^) = (117) 

((M n )^(M„0 M v) = ^2^<U'2Wv'Ai • (118) 

Equations (|115j) to (|1 18j) define an ltd type stochastic process. It is a surprising feature that 
an Ito rather than a Stratonovich type equation [37] arises as a stochastic approximation to a 
system of continuous differential equations. 

Note that we do not claim that the particles are always simply advected by the flow u(r,t) 
when f2 <C 1 and e<l: in the over damped case this is true, but in the underdamped case they 
are not advected. The statement is that the particle separations behave as if the particles were 
being advected, according to the simple model (j!15|) . 



8.2 Lyapunov exponents of Brownian advection model 

The Lyapunov exponents for processes of the type defined by (|115p to (|118p were obtained by 
LeJan |13j . Here we calculate these Lyapunov exponents using a different and much simpler 
approach, facilitating direct comparison with the coefficients in equation ([90]) . 



First Lyapunov exponent 



The Lyapunov exponent Ai is obtained by considering the magnitude of the separations 5r{t) 
\8r(t)\, at successive time steps: 



(119) 



From (I116p we obtain (from here on we drop the subscript labeling the time step on the matrix 
M n ) 

5r 2 (t + At) = 5r 2 (t) n x • (I + M + M T + M T M) m (120) 



where ni is a unit vector such that 5r = 5m\ . Anticipating that two other unit vectors will be 
defined in due course, and writing 

^■Mn„ (121) 



M' 



we obtain 



log c 



Sr(t + At) 
5r(t) 

5r(t + At) 
6r(t) 



1 + 1M' 1X + (M T M)' n 



1 + M' n - \M'\ X + i(M T M)' n + 0(M 3 ) 
M' n - M' 2 n + i(M T M)' n + 0(M 3 ) 
M[ 1 + ^M%-M' 2 11 + 0(M 3 ) . 



(122) 



Taking the expectation value using (JTT7J) and (|118p we obtain 



Ai 



m 2 j 2 



2V 



mi 



(123) 
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where d is the number of space dimensions. Finally we relate this expression for the Lyapunov 
exponent to the results of section [6J showing agreement with the lowest-order terms of the 
expansions e. Using the properties of the elements T>^ Ufl ' u / discussed in section HJ (I123D gives 

A i=^^ i 3(2r-l)=7e 2 (2r-l) (124) 

(here we used the fact that P1313 = P1212 = ^2 = T^i)- This agrees with the leading-order 
term in the expansion (|59"|) . (|9U|) . 



Second Lyapunov exponent 

For the second Lyapunov exponent we follow an approach analogous to that of section [3J we 
consider two vectors 6r± = 5mi and <5r 2 = Sri + <5rc)#n 2 , where n 2 is orthogonal to ni and 
where 59 <C 1. A third unit vector is defined by n3 = ni A nj. The area of the parallelogram 
spanned by the vectors 8r\ and <5r 2 is 6A(t) = \5ri(t) A<5r 2 (t)|. The second Lyapunov exponent 
may be obtained from the relation 



We find 8A(t) = 5r 2 56 and 



5A(t) 

SA(t + At) = |m An 2 + z\ (126) 

where 

z = Mni A n 2 + ni A Mn 2 + Mni A Mn 2 . (127) 

This implies 

It is convenient to express these vectors in terms of components: for example we have Mni = 
M^ni + M 2i n 2 + M^n 3 so that Mni A n 2 = M{in 3 - MLn^ We find n 3 • z = M' n + M' 22 + 
M ii M 22 - A^2i M i2 and z ■ z = M'^ + M'\ 2 + M' 2 n + M' 22 + 2M[ 1 M 2 \ 2 + <3(M 3 ). This gives 



( M ^|t) At) ) 2 = 1 + 2 ( M ii + M 22) + 2M n M 22-2M^M( 2 
+M'l 1 + M% + (M(i + M'22? + 0(M 3 ) 



8A{t + At) 

SA{t) 
5A(t + At) 



l + M' xl + M'22 ~ M' 2l M' l2 + M[ X M^ + \M'\ X + \M'\ 2 + 0(M 3 ) 



log c ^{fy = M n+ M 22-^ 2 M^-lM'?i-lM^ + lMl + iM^ 2 (129) 
+0(M 3 ) . 

Using (|125p . (|117p and (|118p . the sum of the first two Lyapunov exponents is found to be 

Ai + A 2 = — 5— 5-(X>3131 + £>3232 — 2>L111 ~ ^2222 ~ 2£> 2 n 2 ) . (130) 

Expressing (1130j) in terms of the notation used in section [6j we find 

Ai + A 2 = -^(2P 2 - 2Pi - 2P 3 ) = (2r - 2 - 2a) = 7 e 2 (3r - 3) . (131) 

This gives A 2 = 7e 2 (r — 2), which is consistent with the leading-order term of the expansion of 
A 2 in PJ. 
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Third Lyapunov exponent 

The third Lyapunov exponent is determined by considering the sum of the first three Lyapunov 
exponents, which are obtained from the mean logarithmic derivative of the volume element 
SV(t): 

We have 

det(I + M) 

log c det(I + M) = tr log c (I + M) 
tr[M - ±M 2 ] + 0(M 3 ) . (133) 
The sum of the first three Lyapunov exponents is therefore 

Ai + A 2 + A 3 = -i^tr(M 2 )^ 

= r^Pim +V 1221+^1331]- (134) 

In the notation of section [6j (|134p gives 

Ai + A 2 + A 3 = (Z>i + 2P 3 ) = -37£ 2 (1 + 2a) = -3(2 - r) 7 e 2 (135) 

Subtracting (|131 1) we find A3 = — 3 7 e 2 , which agrees with the leading term of (f90|) . 
Equation (|134p may also be written in the form 

1 r°° 

A1 + A2 + A3 = -- / dt (divti>(0,t)divto(O,0)) . (136) 

Here w(r,t) is related to the velocity field w n {r) in equation (|115p by w n (r) = w(r,nAt). A 
similar expression has been quoted by Balkovsky et al. [T7] . 

9 Correction to overdamped limit 

In the previous section, the advective approximation of the solution of equation (J3j) was consid- 
ered. We showed that whenever e <C 1, at leading order in k we can apply a model in which the 
particles are advected by a random flow. The only case where the advective approximation is in- 
adequate in the overdamped limit is the following. Consider the quantity A, defined in equation 
(TT]). Its sign and magnitude determine whether or not clustering occurs, and to which extent. In 
the case of an incompressible flow, A vanishes in the advective approximation, because the flow 
is volume preserving (this also follows from equation (|135p , noting that T = 2 for incompressible 
flow). In order to compute whether A is positive or negative (which determines whether or not 
clustering happens) it is thus necessary to go beyond the advective approximation, as explained 
in the following. 

We consider an incompressible flow field u(r,t). We make use of a result due to Maxey [JJ 
who has derived a 'synthetic' velocity field including a correction term which is proportional 
to I/7 and shown that this modified field describes, to leading order, corrections to the purely 
advective limit. Maxey's synthetic field was subsequently employed in [TBI ESI El to study 
particle-density fluctuations in incompressible fluids in the near-ad vective limit. 



V(t + At) 

W) 

f V(t + At) \ 
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As pointed out by Maxey, this correction term gives rise to a potential component in the 
synthetic advection field which, according to equation (|134p . implies that A > 0. The aim of this 
section is to evaluate this correction, and to express it in terms of the dimensionless parameters 
k and f2. 

Integrating ([3]) one obtains 



Pit) 



I dh exp[-7(t-ti)]/(r(*i),*i) 

J — oo 
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f(Ht),t) 



1 dti eKp[-7(t-ti)]-^-/(r(ti),ti) 

-oo OX\ 



(137) 



ignoring an initial transient. Using / = m^/u this gives, to leading order in Q 1 = (77*) 1 



r = u(r(t),t) 
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u(r(t),t) + (u(r(t),t)-V)r 



(138) 



As pointed out by Maxey (equation (5.9) in pQ) this equation describes advection in a synthetic 
velocity field 



1 

v = u 

7 



- + (..V) W 



(139) 



The sum of the Lyapunov exponents for this advection problem are given by (|136|) . Since 
divu = 7~ 1 V(it • V)u, we have 



Ai + Ao + A3 



(140) 



We will estimate this expression for a particular example, where the correlation function C(R, t) 
in (|35p factorises into a spatial and a time-dependent part: C(R,t) = Co(R)$(t). In the 
incompressible case where V • u = 0, we obtain 



Ai + A2 + A2 



E 



dr u dru 



))w£**®- <141) 



In terms of the characteristic quantities of the flow, we therefore estimate A1+A2+A3 ~ k 4 / (Q 2 t). 
Together with (I124|) this gives, for SI > 1, 



n 2 ' 



(142) 



This result should be contrasted with (|95p . valid for Q <C 1, which implies that A ~ k 2 /£1. Thus 
A is given in terms if different combinations of dimensionless parameters k and f2 for Q <C 1 and 
f2 3> 1. In the underdamped case, A is obtained as a perturbation expansion in e = KfT 1 / 2 . In 
the overdamped case, by contrast, it is obtained as a perturbation expansion in k/Q and k/Q, 2 . 
We note that expression (|142p is always small in the domain of its validity: £1 3> 1 and k < 1. 
The results of section [6] by contrast, remain valid when A is of order unity. 



10 Concluding remarks 
10.1 Conclusions 

In this paper we have characterised the local properties of the motion of inertial particles in 
a short-time correlated random flow model. We calculated the Lyapunov exponents using an 
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exact series expansion, and characterised the rate of caustic formation in terms of an escape 
process. 

The results give important new insights into the mechanism for particle clustering. The 
relevant dimensionless parameter is e ~ kJI -1 / 2 , and we find that the dimension deficit is 
significant when e is of order unity, reaching a maximum value of A max 0.35 at e max as 
0.21. This contradicts the accepted explanation for particle clustering, based on the centrifuge 
mechanism, because (when n <C 1) it implies that clustering onto a fractal set can occur when 
the Stokes number O -1 is large. For fully developed turbulent flow, we have k = O(l) so that 
e ~ SI -1 / 2 , and in that case our theory does predict that clustering onto a fractal set occurs 
when the Stokes number is of order unity. This is consistent with numerical experiments. We 
cannot exclude the possibility that the centrifuge effect has some relevance, however we note 
that our theory predicts A max = 0.35 for the maximal dimension deficit, see figure [TJ}, also [5J. 
The value obtained from direct numerical simulations of a Navier-Stokes flow is very similar, 
A max = 0.4 |]. 

For e > e c ~ 0.33, the dimension deficit is zero, and the particles will not cluster onto a 
fractal set. However, as e increases the rate of caustic formation increases. The particle density 
can diverge on caustic lines, as described in |20j. 

Figure 7 is a schematic illustration of how the unmixing of particles suspended in an incom- 
pressible fluid (that is a flow with T = 2) depends upon the parameters k and £1. Instead of 
clustering being confined to $7 ~ 1, it may occur on a ray in logarithmic parameter space ex- 
tending into the underdamped region. The solid line in figure [HJ e ~ kQ^ 1 / 2 = const., indicates 
schematically the critical line where A reaches zero. Above the solid line A is always positive, 
but tends to zero for small e as A = 10e 2 ~ k 2 /Q. In section 8 we showed that in the limit of 
e — > the dynamics behaves like advection, even if Q -C 1. 

10.2 Experimental considerations 

We make some final comments concerning the situations in which clustering should be observable. 
The velocity field of fully-developed turbulent flow has a power-law energy spectrum, with an 
upper and lower cutoff [26]. The smaller length scale is the Kolmogorov length, defined by 
r\ = (v 3 /E) l / A where v is the kinematic viscosity and £ is the rate of dissipation per unit mass 
of fluid. This is the size of the smallest vortices generated by the turbulence. In section 14.31 we 
show that it is this small length scale r] which corresponds to the correlation length £ in our 
theory. Our results indicate that clustering can certainly occur on a scale smaller than £, but 
what happens on larger length scales is less certain. 

Lengthscale for clustering 

It is not entirely clear whether particles in turbulent flows can exhibit clustering effects on 
length scales greater than the Kolmogorov length, rj. Particles separated by distances greater 
than rj will be separated by Richardson diffusion [39], but that is not necessarily incompatible 
with some type of clustering. Numerical evidence for such an effect has been presented ]7J, but 
the range of Stokes and Reynolds numbers are not sufficiently large to allow clear conclusions. 
At very small Stokes number, suspended particles are advected and must be mixed to uniform 
concentration in a turbulent flow. At large Stokes number, however, it is conceivable that 
particles could be centrifuged away from vortices with length scale L 3> r/ when their rotation 
period ti is such that the length-dependent Stokes number St/, = 1/(^)T£) is of order unity [40] . 
The Kolmogorov theory gives ~ (L 2 /^) 1 / 3 . Consequently this picture implies clustering on a 
length scale L ~ (£ /7 3 ) 1 / 2 . 

However, there is an argument which suggests that such large-scale clusters might not form. 
At large Stokes number, the formation of numerous caustics implies that the velocity field of 
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the suspended particles is multi-valued, so that in the limit St — > oo the suspended particles are 
expected to behave as a gas [H]. Particles at the same position move with relative speed Av; 
using the Kolmogorov cascade principle, we estimate that the variance of Av is (Av 2 ) ~ 6/j 
[42] . The motion of the suspended particles is damped at a rate 7, so that they travel in an 
approximately straight line for a distance Av/j. We see that Af/7 ~ (1S/7 2 ) 1 / 2 ~ L. So 
particles within a cluster of size L travel distances of order L in random directions before their 
relative motion is damped out. This argument suggests that clustering may be difficult to 
observe on scales larger than the Kolmogorov length. 



Effects of gravity 

The Kolmogorov scale is rather insensitive to the value of £: for v ~ 10~ 5 m 2 s~ 1 (air at standard 
atmospheric conditions) and £ ~ lm 2 s -3 , we have r\ ~ 10~ 4 m. Our results concerning the 
Lyapunov exponents are relevant to clustering of particles on length scales below the Kolmogorov 
length and therefore describe clustering on very small length scales. They would be relevant to 
the aggregation of aerosol particles, but it seems unlikely that it could explain the macroscopic 
clustering observed in the experiments of Fessler and Eaton and [2] . 

In discussing the conditions under which the clustering effect is observable, we must also 
consider the effects of gravity. Gravitational effects may be assumed negligible when motion 
at the terminal velocity moves particles by a distance which is small compared to £ in time r. 
Noting that the terminal velocity is g/7, this condition may be written gr 2 /£f2 <C 1. For fully 
developed turbulence, using the Kolmogorov theory we estimate 

This quantity is small only for highly energetic turbulence. Thus, in terrestrial applications, 
direct application of the model to three-dimensional turbulent flows requires very intense tur- 
bulence if the effects of gravity are to be neglected. 
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Figure 1: a Dimension deficit A as a function of e, in the limit k — > 0; data from simulation 
of our short-time correlated random-flow model obtained as described in section [3721 (o). from a 
simulation of the Langevin equation (1441 - 148ft . solid red line, and results of a Borel summation 

of our perturbation series, using a Pade approximant of order P 3 / 3 ( ). The figure shows 

that there is pronounced clustering when e = 0(1): in the limit as k — > this happens when 
St 3> 1. b Dimension deficit A as a function of Stokes number St; results (taken from [3]) of 
direct numerical simulations of particles embedded in Navier-Stokes flows (□), and results from 
Langevin equations of our short-time correlated random-flow model, (|44l - H8]). solid red line. 
The latter are the same as those in figure [1^, with St = e 2 / k 2 and k = 0.256 was chosen to give 
a good fit between the curves, as judged by the eye. 



X 




Figure 2: Schematic illustration of a caustic in one spatial dimension: the velocity x of the 
particles as a function of position x is initially single- valued (curve (a)) but particles with a 
large velocity overtake slower-moving particles, so that at a later time the particle velocity is 
multi- valued (curve (b)). The region where the velocity is a multi- valued function is bounded 
by two fold caustics. 
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Figure 3: Lyapunov exponents Xj as a function of e for j = 1,2,3 and three different values of T. 
Panels a - c show numerical simulations using a discretised form of the equations of motion ([3]), 
replacing the force field by a random variable with appropriate statistics as described in section 
I3.2t o (j = 1), □ (j = 2), and o (j = 3). Also shown are results of Monte-Carlo simulation of 
the Langevin equation (f44j) . solid lines, as well as the lowest-order estimates (f94j) . dashed lines. 
Panel a: T = 2 (incompressible flow). Panel b: T = 1/3 (pure potential flow). Panel c: T = 1. 
Panel d shows large-e behaviour of the numerical data from a, also shown are fits to the large-e 
behaviour (|49p . dashed lines. 
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Figure 4: Illustrating the growth of the perturbation coefficients ch with increasing order n for 
the first Lyapunov exponent, T = 2. Shown is the ratio c^+i/oh^ as a function of n. The solid 

line is the asymptotic relation c£L\/dn ~~ " ' " J 
constant C. 



6n, consistent with Cn ~ C 6 n (n — 1)!, for some 
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Figure 5: Borel-Pade summation of the perturbation expansion of Xj. a: Numerical data for Ai 
as a function of e for T = 2 (o) (obtained as described in section 13, 2 j) . compared with results 

of Borel-Pade summation, equations ([92]) and ([93]) ; P 5 / 5 ( ), fWio (— • — • — ), P15/15 

(_.._.. _) f anc [ P 2 q/2o ( )• and -P23/23 ( )• b: Numerical data for Ai (o) and A2 

(□) as a function of e for T = 1. Also shown are results of Pade-Borel summation (using a Pade 
approximant -P23/23) f° r J = 1 (solid line) and j = 2 (dashed line). Since according to equation 

(pi) c} 2) = -c{ 1} at T = 1, those curves are identical apart from an overall minus sign. But 
using the relation Ai + A2 ~ Cexp[— 1/(6 e 2 )], with Ai obtained by Borel-Pade summation and 
C = 0.79 a fitted parameter (—• — •—) gives a much better fit to A2 than the dashed line. 
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Figure 6: a dependence of the Lyapunov exponents on e, for T = 2: the symbols show numerical 
data obtained as described in section 13.21 and the solid lines are the results of Borel-Pade 
summation. The results for A3 depend upon the number of orders of the perturbation series 



included in the Borel sum, P, 
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and P 
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dimension deficit obtained from the curves in figure [6k using equation ([T]): the same line styles 
are used for different Pade approximants. The solid red curve shows the results of Monte-Carlo 
simulations of the Langevin equation ([441 - l48l) , and is the same as in figure [Th. The dashed red 
line is the asymptotic approximation (|95p . 
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Figure 7: Rate of caustic formation J' = J/7 as a function of e, for T = 2, Monte-Carlo 
simulations of the Langevin equation (|44l - EE8|) . solid red line, and numerical data obtained as 
described in section [331 (o). The results are asymptotic to J' = C exp(— $ caust i c /e 2 ) in the limit 
as e — > 0, with C ~ 0.28 and ^caustic ~ 0.125 (dashed line). 




Figure 8: Schematic diagram showing how the behaviour of particles suspended in an incom- 
pressible fluid depends upon the parameters k and 0. 
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